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In this paper we present the results ol A-body simulations with a scalar field coupled differently 
to cold dark matter (CDM) and baryons. The scalar field potential and coupling function are chosen 
such that the scalar field acquires a heavy mass in regions with high CDM density and thus behaves 
like a chameleon. We focus on how the existence of the scalar field affects the formation of nonlinear 
large-scale structure, and how the different couplings of the scalar field to baryons and CDM particles 
lead to different distributions and evolutions for these two matter species, both on large scales and 
inside virialized halos. As expected, the baryon-CDM segregation increases in regions where the fifth 
force is strong, and little segregation in dense regions. We also introduce an approximation method 
to identify the virialized halos in coupled scalar field models which takes into account the scalar 
field coupling and which is easy to implement numerically. It is find that the chameleon nature of 
the scalar field makes the internal density profiles of halos dependent on the environment in a very 
nontrivial way. 

PACS numbers: 04.50.Kd 



I. INTRODUCTION 

The origin and nature of the dark energy [1] is one of 
the most difficult challenges facing physicists and cosmol- 
ogists now. Among all the proposed models to tackle this 
problem, a scalar held is perhaps the most popular one 
up to now. The scalar held, denoted by (p, might only in- 
teract with other matter species through gravity, or have 
a coupling to normal matter and therefore producing a 
fifth force on matter particles. This latter idea has seen 
a lot of interests in recent years, in the light that such a 
coupling could potentially alleviate the coincidence prob- 
lem of dark energy [2] and that it is commonly predicted 
by low energy effective theories from a fundamental the- 
ory. 

Nevertheless, if there is a coupling between the scalar 
held and baryonic particles, then stringent experimental 
constraints might be placed on the fifth force on the lat- 
ter provided that the scalar held mass is very light (which 
is needed for the dark energy). Such constraints severely 
limit the viable parameter space of the model. Different 
ways out of the problem have been proposed, of which the 
simplest one is to have the scalar field coupling to dark 
matter only but not to standard model particles, there- 
fore evading those constraints entirely. This is certainly 
possible, especially because both dark matter and dark 
energy are unknown to us and they may well have a com- 
mon origin. Another interesting possibility is to have the 
chameleon mechanism [3-6] , by virtue of which the scalar 
field acquires a large mass in high density regions and 
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thus the fifth force becomes undetectablly short-ranged, 
and so also evades the constraints. 

Study of the cosmological effect of a chameleon scalar 
field shows that the fifth force is so short-ranged that it 
has negligible effect in the large scale structure formation 
[7] for certain choices of the scalar field potential. But it 
is possible that the scalar held has a large enough mass in 
the solar system to pass any constraints, and at the same 
time has a low enough mass (thus long range forces) on 
cosmological scales, producing interesting phenomenon 
in the structure formation. This is the case of some f(R) 
gravity models [8, 9], which survives solar system tests 
thanks again to the chameleon effect [10-13]. Note that 
the f{R) gravity model is mathematically equivalent to 
a scalar field model with matter coupling. 

No matter whether the scalar field couples with dark 
matter only or with all matter species, it is of general in- 
terests to study its effects in cosmology, especially in the 
large scale structure formation. Indeed, at the linear per- 
turbation level there have been a lot of studies about the 
coupled scalar held and f(R) gravity models which en- 
able us to have a much clearer picture about their behav- 
iors now. But linear perturbation studies do not conclude 
the whole story, because it is well known that the matter 
distribution at late times becomes nonlinear, making the 
behavior of the scalar field more complex and the linear 
analysis insufficient to produce accurate results to con- 
front with observations. For the latter purpose the best 
way is to perform full A^-body simulations [14] to evolve 
the individual particles step by step. 

A^-body simulations for scalar held and relevant mod- 
els have been performed before [15-22]. For example, in 
[21] the simulation is about a specific coupled scalar field 
model. This study however does not obtain a full solu- 
tion to the spatial conhguration of the scalar field, but 
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instead simplifies the simulation by assuming that the 
scalar field's effect is to change the value of the gravi- 
tational constant, and presenting an justifying argument 
for such an approximation. As discussed in [23, 24], this 
approximation is only good in certain parameter spaces 
and for certain choices of the scalar field potential, and 
therefore full simulations tare needed to study the scalar 
held behaviour more rigorously. 

Recently there have also appeared A-body simulations 
of the f(R) gravity model [25, 26], which do solve the 
scalar degree of freedom explicitly. However, embedded 
in the f(R) framework there are some limitations in the 
generality of these works. As a first thing, f(R) gravity 
model (no matter what the form / is) only corresponds 
to the couple scalar field models for a specific value of 
coupling strength [27]. Second, in f(R) models the cor- 
rection to standard general relativity is through the mod- 
ification to the Poisson equation and thus to the gravi- 
tational potential as a whole [25], while in the coupled 
scalar field models we could clearly separate the scalar 
fifth force from gravity and analyze the former directly 
[23]. Also, in f(R) models, as well as the scalar-tensor 
theories, the coupling between matter and the scalar field 
is universal (the same to dark matter and baryons), while 
in the couple scalar field models it is straightforward to 
switch on/off the coupling to baryons and study the ef- 
fects on baryonic and dark matter clusterings respectively 
(as we will do in this paper). Correspondingly, the gen- 
eral framework of A-body simulations in coupled scalar 
field models could also handle the situation where the 
chameleon effect is absent and/or scalar field only couples 
to dark matter, and thus provide a testbed for possible 
violations of the weak equivalence principle. 

In this paper we shall go beyond [23] and consider the 
case where the chameleon scalar field couples differently 
to different species of matter. To be explicit, we consider 
two matter species, and let one of them have no coupling 
to the scalar field. Because it is commonly believed that 
normal baryons, being observable in a variety of experi- 
ments, should have extremely weak (if any) coupling to 
scalar fields, we call the uncoupled matter species in our 
simulation " baryons" . It is however reminded here that 
this matter species is not really baryonic in the sense 
that it does not experience normal baryonic interactions. 
The inclusion of true baryons will make the investigation 
more complicated and is thus beyond the scope of the 
present work. 

The paper is organized as follows: in § II we list the 
essential equations to be implemented in the A-body sim- 
ulations and describe briefly the difference from normal 
LCDM simulations. § III is the main body of the paper, 
in which § III A gives the details about our simulations, 
such as code description and parameter set-up, § IIIB 
displays some preliminary results for visualization, such 
as baryon/CDM distribution, potential/scalar field con- 
figuration and the correlation between the fifth force (for 
CDM particles) and gravity; § III C quantifies the nonlin- 
ear matter power spectrum of our model, especially the 



difference from LCDM results and the bias between CDM 
and baryons; § III D briefly describes the essential mod- 
ifications one must bear in mind when identifying viri- 
lized halos from the simulation outputs and shows the 
mass functions for our models; § III E we pick out two 
halos from our simulation box and analyzes their total 
internal profiles, as well as their baryonic/CDM density 
profiles. We finally summarize in § IV. 



II. THE EQUATIONS 

In this section we first describe the method to simulate 
structure formation with two differently coupled matter 
species and the appropriate equations to be used. Those 
equations for a single matter species have been discussed 
in details previously in [23] , but the inclusion of different 
matter species requires further modifications and we list 
all these for completeness. 



A. The Fundamental Equations 



The Lagrangian for our coupled scalar field model is 



C = 



1 



R 



- V a <pV a <p 



+ V( V )-C(<p)Lcvu + Cs(l) 



where R is the Ricci scalar, k = 8ttG with G Newton's 
constant, <p is the scalar field, V(tp) is its potential energy 
and C(y>) its coupling to dark matter, which is assumed 
to be cold and described by the Lagrangian £cdm ■ £s in- 
cludes all other matter species, in particular our baryons. 
The contribution from photons and neutrinos in the A- 
body simulations (for late times, i.e., z ~ 0(1)) is negli- 
gible, but should be included when generating the matter 
power spectrum from which the initial conditions for our 
A-body simulations are obtained (see below). 

The dark matter Lagrangian for a point-like particle 
with bare mass mo is 



-CcDM(y) 



m 



■Ay - xo)\Jgabi%x b (2) 



where y is the coordinate and Xo is the coordinate of the 
centre of the particle. From this equation it can be easily 
derived that 



rpab 
J CDM 



m 



5{y-x )x%x b 



(3) 



Also, because g a bi C Qx\ = g a t,u a u = 1 where u a is the 
four velocity of the dark matter particle, the Lagrangian 
could be rewritten as 



•CcDM(y) — 
which will be used below. 



m 



(4) 
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Eq. (3) is just the energy momentum tensor for a single 
dark matter particle. For a fluid with many particles the 
energy momentum tensor will be 



rpab 
1 CDM 



PCDMU a U b , 



(5) 



in which V is a volume microscopically large and macro- 
scopically small, and we have extended the 3-dimcnsional 
d function to a 4-dimensional one by adding a time com- 
ponent. Here u a is the averaged 4- velocity of the collec- 
tion of particles inside this volume, and is not necessarily 
the same as the 4- velocity of the observer. 
Meanwhile, using 



rpab 



Sg a b 



(6) 



it is straightforward to show that the energy momentum 
tensor for the scalar field is given by 



-^ipab 



= wvV-5 



(7) 



So the total energy momentum tensor is 



Tab 



Va^Vf,^ - g a b 



1 ab 



(8) 



where T^, DM = PcnMU a Ub, T^ b is the energy momentum 
tensor for all other matter species including baryons, and 
the Einstein equation is 



Gab 



nT a b 



(9) 



where G a b is the Einstein tensor. Note that due to the 
coupling between the scalar field ip and the dark mat- 
ter, the energy momentum tensors for either will not be 
conserved, and we have 



CDMab 



C Af) (ab r 

■ C{<p) [g 



CDM 



T 



CDMab\ 



Vb^lO) 



where throughout this paper wc shall use a v to denote 
the derivative with respect to ip. 

Finally, the scalar field equation of motion (EOM) from 
the given Lagrangian is 



dV(ip) dC{ip) 



£CDM 



dip dip 
where □ = V a V Q . Using Eq. (4) it can be rewritten as 



n<p- 



9V(<P) 
dip 



PCDM" 



dip 



0. 



(11) 



Eqs. (8, 9, 10, 11) summarize all the physics that will 
be used in our analysis. 



We will consider a special form for the scalar field po- 
tential, 



V(ip) 



V 



[1 - exp (/3^ip)} 



(12) 



where /i and /? are dimensionless constants while Vq has 
mass dimension four. As has been discussed in [23], 
to evade observational constraints and /? can be set to — 1 
without loss of generality, since we can always rescale ip 
as we wish. Meanwhile, the coupling between the scalar 
field and dark matter particle is chosen as 



C(ip) = cxp(jy/Kip), 



(13) 



where 7 > is yet another dimensionless constant char- 
acterizing the strength of the coupling. 

As discussed in [23] , the two dimensionless parameters 
(i and 7 have clear physical meanings: roughly speaking, 
jj, controls the time when the scalar field becomes im- 
portant in cosmology while 7 determines how important 
the scalar field would ultimately be. In fact, the potential 
given in Eq. (12) is partly motivated by the f(R) cosmol- 
ogy [11], in which the extra degree of freedom behaves as 
a coupled scalar field in the Einstein frame. As we can 
see from Eq. (12), the potential V — , 00 when ip — »• 
while V —> Vq when ip — > 00. In the latter case, however, 
C — >• 00, so that the effective total potential 



V efj(<P) = V(<p) + PCBMC(ip) 



(14) 



has a global minimum at some finite ip. If the total poten- 
tial V e ff(ip) is steep enough around this minimum, then 
the scalar field becomes very heavy and thus follows its 
minimum dynamically, as is in the case of the chameleon 
cosmology (see e.g. [7]). If V e ff is not steep enough at 
the minimum, however, the scalar field will experience 
a more complicated evolution. These two different cases 
can be obtained by choosing appropriate values of 7 and 
fi: if 7 is very large or \x is small then we run into the for- 
mer situation and if 7 is small and /1 is large we have the 
second. In reality, the situation can get even more com- 
plicated because when 7, which characterizes the cou- 
pling strength, increases, the CDM evolution could also 
get severely affected, which in turn has back-reactions on 
the scalar field itself. 



B. Nonrelativistic Limit 

The iV-body simulation only probes the motion of par- 
ticles at late times, and we are not interested in extreme 
conditions such as black hole formation/evolution, which 
mean that taking the non-relativistic limit of the above 
equations should be a sufficient approximation for our 
purpose. 

The existence of the scalar field and its (different) cou- 
plings to matter particles lead to the following changes 
to the ACDM model: Firstly, the energy momentum ten- 
sor has a new piece of contribution from the scalar field; 
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secondly, the energy density of dark matter in gravita- 
tional field equations is multiplied by the function C(tp), 
which is because the coupling to scalar field essentially 
renormalizes the mass of dark matter particles; thirdly, 
dark matter particles will not follow geodesies in their 
motions as in ACDM, but rather the total force on them 
has a contribution, the fifth force, from the exchange of 
scalar field quanta; finally, CDM particles must be dis- 
tinguished from baryons so that the fifth force only acts 
on the former and these two species only interact gravi- 
tationally. This last point is one main difference between 
the present work and a previous one [23]. 

These imply that the following things need to be mod- 
ified or added: 

1. The scalar field tp equation of motion, which deter- 
mines the value of the scalar field at any given time 
and position; 

2. The Poisson equation, which determines the grav- 
itational potential (and thus gravity) at any given 
time and position, according to the local energy 
density and pressure, which include the contribu- 
tion from the scalar field (as obtained from tp equa- 
tion of motion); 

3. The total force on the dark matter particles, which 
is determined by the spatial configuration of tp, just 
like gravity is determined by the spatial configura- 
tion of the gravitational potential; 

4. The CDM and baryonic particles must be tagged 
respectively so that the code knows to assign forces 
correctly to different species. 

We shall describe these one by one now. 

For the scalar field equation of motion, we denote tp as 
the background value of tp and Sip = tp — if as the scalar 
field perturbation. Then Eq. (11) could be rewritten as 



S'tp + 3HStp + V 2 r ip + V lV (tp) - V tV (<p) 

+PCDmC iV (<^) - p~CDMC }tp (tf) 







by subtracting the corresponding background equation 
from it. Here V ra is the covariant spatial derivative with 
respect to the physical coordinate r = ax with x the con- 
formal coordinate, and V 2 = V ra V°. V ra is essentially 
the V a , but because here we are working in the weak field 
limit we approximate it as = — (d 2 ^ + d% y + d' 2 J by 

assuming a flat background; the minus sign is because our 
metric convention is (+, — , — , — ) instead of (— , +, +, +). 
For the simulation here we will also work in the quasi- 
static limit, assuming that the spatial gradient is much 
larger than the time derivative, |V r <y9| ^> |^| (which will 
be justified below). Thus the above equation can be fur- 
ther simplified as 



- 2 dl(a5tp) 



in which d 2 = — V 2 = + (d 2 + d 2 + d 2 ) is with respect 
to the conformal coordinate x so that V x = aV r , and 
we have restored the factor c 2 in front of V 2 (the tp here 
and in the remaining of this paper is c~ 2 times the tp 
in the original Lagrangian unless otherwise stated). Note 
that here V and /?cdm both have the dimension of mass 
density rather than energy density. 

Next look at the Poisson equation, which is obtained 
from the Einstein equation in weak-field and slow-motion 
limits. Here the metric could be written as 

ds 2 = (1 + 2<j))dt 2 - (1 - 2i>)S. lJ dr l dr : > (16) 

from which we find that the time-time component of the 
Ricci curvature tensor R° = — V 2 t/>, and then the Ein- 
stein equation R a t = k (T a b — \g a bT^j gives 



R° = -v 2 < 



= 2 (^tot + 3 Ptot) (17) 



where ptot and j>tot are respectively the total energy 
density and pressure. The quantity V 2 </> can be expressed 
in terms of the comoving coordinate x as 

1 ^.o 1 



(18) 



V> = — V£ - - -aax 
a z \a 2 

= ^V 2 <&-3^ 
a 3 x a 

where we have defined a new Newtonian potential 

$ = a<p + ia 2 ax 2 
and used Vjx 2 = 6. Thus 



(19) 



(20) 



2 (ptot + 3ptot) — 2 (/'tot + 3ptot) 



where in the second step we have used Eq. (17) and the 
Raychaudhrui equation, and an ovcrbar means the back- 
ground value of a quantity. Because the energy momen- 
tum tensor for the scalar field is given by Eq. (7), it is 
easy to show that p v + 3p v = 2 \tp 2 — V{tp)\ and so 

V 2 $ 

= -4irGa 3 {pcbuC(p) +p B + 2[p 2 - V(tp)] } 
+4nGa 3 {p CDM C(^) + p B + 2 [<f 2 - V{tp)} } . 



Now in this equation ip 2 — tp 2 = 2tp8tp + 5p <C (V r p) 2 
in the quasi-static limit and so could be dropped safely. 
So we finally have 



° 3 [ v ,<p(f) ~ V^iv) + PcbmC iV {p) - pcbmC iV (p)} 



(15) d 2 $ = 47rGa 3 [pcdmC^) - pcdmC(^)] 



+4irGa 3 [p B - p B ] - 8^Ga 3 [V(tp) - V(<p)]{21) 
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Finally, for the equation of motion of the dark matter 
particle, consider Eq. (10). Using Eqs. (3, 4), this can be 
reduced to 



• r h-'V 



ab 



VbP- (22) 



Obviously the left hand side is the conventional geodesic 
equation and the right hand side is the new fifth force 
due to the coupling to the scalar field. Note that be- 
cause g ab — u a u b = h ab is the projection tensor that 
projects any 4-tensor into the 3-space perpendicular to 
u a , so (g ab — « a u b ) V a = V b is the spatial derivative in 
the 3-space of the observer and perpendicular to u a \ con- 



sequently the fifth force 



cy(y) 

C(<p) 



V a </3 = V a logC(ip) has no 



component parallel to u a (the time component), indicat- 
ing that the energy density of CDM will be conserved and 
only the particle trajectories are modified, as mentioned 
in [23]. Remember that u a in Eq. (22) is the 4- velocity 
of individual particles, but from Eq. (15) we see that Sip 
is computed in the fundamental observer's frame (where 
density perturbation is calculated), so if we also want to 
work on Eq. (22) in the fundamental observer's frame (so 
that we can use the Sip from Eq. (15) directly), then we 
must rewrite Eq. (22) by substituting 

(g ab - u a u b ) V b ip = (g ab -u a u b - u a v b -u b v a )V b ip 
« (g ab - ii a u b ) V b <p - <pv a 

up to first order in perturbations, in which u a is the 4- 
velocity of the fundamental observer and v is the peculiar 
velocity of the particle. Then the first term in the above 
expression is the gradient of Sip observed by the funda- 
mental observer (rather than an observer comoving with 
the particle) and the second term is a velocity dependent 
acceleration [21]. In [22] it is claimed that the second 
term is of big importance; in our simulations, however, 
this term will be neglected (from here on) because it de- 
pends on ip, which is very small due to the chameleon 
nature of the model. We have checked in a linear pertur- 
bation computation that removing this term only changes 
the matter power spectrum by less than 0.0001%. 

Now in the non-relativistic limit the spatial compo- 
nents of Eq. (22) can be written as 



d 2 r 

IT 2 = 



C(p) 



V r ip 



(23) 



where t is the physical time coordinate. If we instead use 
the comoving coordinate x, then this becomes 



a 1 - 
x + 2-x = T V X $- 

a a J 



1 C v (<p)~ 
a 2 C(ip) 



V x¥ ? (24) 



where we have used Eq. (19). The canonical momentum 



conjugate to x is p = a 2 x so we have now 



ax 


P 


~dt 


a 2 ' 


rfPCDM 

dt 


--V x $ 

a 


dp B 

dt 


--v x $, 

a 



C(ip) 



(25) 
(26) 
(27) 



in which Eq. (2G) is for CDM particle and Eq. (27) is for 
baryons. Note that according to Eq. (26) the quantity 
a\og[C{ip)] acts as a new piece of potential: the potential 
for the fifth force. This is an important observation and 
we will come back to it later when we calculate the escape 
velocity of CDM particles within a virializcd halo. 

Eqs. (15, 21, 25, 26, 27) will be used in the code to eval- 
uate the forces on the dark matter particles and evolve 
their positions and momenta in time. 



C. Internal Units 

In our numerical simulation we use a modified version 
of MLAPM ([32], see III A), and we will have to change 
our above equations in accordance with the internal units 
used in that code. Here we briefly summarize the main 
features. 

MLAPM code uses the following internal units (with 
subscript c ): 



x c 

Pc 

t c 
$c 
Pc 



x/B, 

P/(H B) 

tH 

P/P, 



(28) 



in which B is the present size of the simulation box and 
Hq is the present Hubble constant, and p. with subscript, 
can represent the density for either CDM (p c ,Cdm) or 
baryons (p c ,b)- Using these newly-defined quantities, it is 
easy to check that Eqs. (25, 26, 21, 15) could be rewritten 
as 



dx c 


Pc 


dt c 


a z 


dp c 


1 


dt c 


a 


V 2 $ c - 


¥ 




3 




+ 2 



c, 
~c 



c^Vtp 



CDmC I Pc,CDM 



c 



Ob (Pc.b - 1) - k 



V-V 

H 



(29) 
(30) 

(31) 



and 



(BH ) 



■V 2 (aip) 



= -Or 



[ -c, 



- 1 



v v -v 



f „3, 



Hi 



: (32) 
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where ficDM is the present CDM fractional energy den- 
sity, we have again restored the factor c 2 and again the ip 
is c -2 times the ip in the original Lagrangian. Note that 
in Eq. (30) the term in the bracket on the right hand side 
only applies to CDM but not to baryons. Also note that 
from here on we shall use V = d Xc , V 2 = <9 Xr ■ d Xc unless 
otherwise stated, for simplicity. 
We also define 

X = \/k<p, (33) 
u = ln(e x - 1) (34) 



to be used below. 

Making discretized version of the above equations for 
iV-body simulations is non-trivial task. For example, the 
use of variable u instead of ip (Appendix A) helps to pre- 
vent ip < 0, which is unphysical, but numerically possible 
due to discretization. We refer the interested readers to 
Appendix A to the whole treatment, with which we can 
now proceed to do A-body runs. 

III. SIMULATION AND RESULTS 
A. Simulation Details 

1. The N-body Code: MLAPM 

The full name of MLAPM is Multi-Level Adaptive Par- 
ticle Mesh code. As the name has suggested, this code 
uses multilevel grids [33-35] to accelerate the convergence 
of the (nonlinear) Gauss-Seidel relaxation method [34] in 
solving boundary value partial differential equations. But 
more than this, the code is also adaptive, always refin- 
ing the grid in regions where the mass/particle density 
exceeds a certain threshold. Each refinement level form 
a finer grid which the particles will be then (re)linked 
onto and where the field equations will be solved (with a 
smaller time step). Thus MLAPM has two kinds of grids: 
the domain grid which is fixed at the beginning of a sim- 
ulation, and refined grids which are generated according 
to the particle distribution and which are destroyed after 
a complete time step. 

One benefit of such a setup is that in low density re- 
gions where the resolution requirement is not high, less 
time steps arc needed, while the majority of computing 
sources could be used in those few high density regions 
where high resolution is needed to ensure precision. 

Some technical issues must be taken care of however. 
For example, once a refined grid is created, the particles 
in that region will be linked onto it and densities on it 
are calculated, then the coarse-grid values of the grav- 
itational potential are interpolated to obtain the corre- 
sponding values on the finer grid. When the Gauss-Scidcl 
iteration is performed on refined grids, the gravitational 
potential on the boundary nodes are kept constant and 



only those on the interior nodes are updated according to 
Eq. (All): just to ensure consistency between coarse and 
refined grids. This point is also important in the scalar 
field simulation because, like the gravitational potential, 
the scalar field value is also evaluated on and communi- 
cated between multi-grids (note in particular that differ- 
ent boundary conditions lead to different solutions to the 
scalar field equation of motion) . 

In our simulation the domain grid (the finest grid that 
is not a refined grid) has 128 3 nodes, and there are a 
ladder of coarser grids with 64 3 , 32 3 , 16 3 , 8 3 , 4 3 nodes 
respectively. These grids are used for the multi-grid ac- 
celeration of convergence: for the Gauss-Seidel relaxation 
method, the convergence rate is high upon the first sev- 
eral iterations, but quickly becomes very slow then; this 
is because the convergence is only efficient for the high 
frequency (short-range) Fourier modes, while for low fre- 
quency (long-range) modes more iterations just do not 
help much. To accelerate the solution process, one then 
switches to the next coarser grid for which the low fre- 
quency modes of the finer grid are actually high frequency 
ones and thus converge fast. The MLAPM solver adopts 
the self-adaptive scheme: if convergence is achieved on a 
grid, then interpolate the relevant quantities back to the 
finer grid (provided that the latter is not on the refine- 
ments) and solve the equation there again; if convergence 
becomes slow on a grid, then go to the next coarser grid. 
This way it goes indefinitely except when converged solu- 
tion on the domain grid is obtained or when one arrives at 
the coarsest grid (normally with 2 3 nodes) on which the 
equations can be solved exactly using other techniques. 
For our scalar field model, the equations are difficult to 
solve anyway, and so we truncate the coarser-grid series 
at the 4 3 -node one, on which wc simply iterate until con- 
vergence is achieved. Furthermore, we find that with the 
self-adaptive scheme in certain regimes the nonlinear GS 
solver tends to fall into oscillations between coarser and 
finer grids; to avoid such situations, we then use V-cycle 
[34] instead. 

For the refined grids the method is different: here one 
just iterate Eq. (All) until convergence, without resort- 
ing to coarser grids for acceleration. 

As is normal in the Gauss-Scidcl relaxation method, 
convergence is deemed to be achieved when the numerical 
solution after n iterations on grid k satisfies that the 
norm || • || (mean or maximum value on a grid) of the 
residual 

e k = L k (u k )-f k , (36) 

is smaller than the norm of the truncation error 

r k = L k -\Kui)-K[L k {u k n )] (37) 

by a certain amount, or, in the V-cycle case, the reduc- 
tion of residual after a full cycle becomes smaller than a 
predefined threshold (indeed the former is satisfied when- 
ever the latter is). Note here L k is the discretization of 
the differential operator Eq. (A9) on grid k and L k ~ 1 
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a similar discretization on grid k — 1, fk is the source 
term, 1Z is the restriction operator to interpolate values 
from the grid k to the grid k — 1. In the modified code 
we have used the full-weighting restriction for 1Z. Corre- 
spondingly there is a prolongation operator V to obtain 
values from grid k — 1 to grid fc, and we use a bilinear 
interpolation for it. For more details see [32]. 

MLAPM calculates the gravitational forces on parti- 
cles by centered difference of the potential 4> and propa- 
gate the forces to locations of particles by the so-called 
triangular-shaped-cloud (TSC) scheme to ensure momen- 
tum conservation on all grids. The TSC scheme is also 
used in the density assignment given the particle distri- 
bution. 

The main modifications to the MLAPM code for our 
model are: 

1. We have added a parallel solver for the scalar field 
based on Eq. (A3). The solver uses a nonlinear 
Gauss-Seidel method and the same criterion for 
convergence as the (linear) Gauss-Seidel Poisson 
solver. 

2. The solved value of u is then used to calculate lo- 
cal mass density and thus the source term for the 
Poisson equation, which is solved using fast Fourier 
transform. 

3. The fifth force is obtained by differentiating the u 
just like the calculation of gravity. 

4. The momenta and positions of particles are then 
updated taking in account of both gravity and the 
fifth force. 

There are a lot of additions and modifications to ensure 
smooth interface and the newly added data structures. 
For the output, as there are multilevel grids all of which 
host particles, the composite grid is inhomogeneous and 
thus we choose to output the positions, momenta of the 
particles, plus the gravity, fifth force and scalar field value 
at the positions of these particles. We can of course easily 
read these data into the code, calculate the corresponding 
quantities on each grid and output them if needed. 

2. Inclusion of Baryons 

As mentioned above, the most important difference of 
the present work from [23] is the inclusion of baryons - 
the particles which do not couple to the scalar field. The 
baryons do not contribute to the scalar field equation of 
motion and are not affected by the scalar fifth force, at 
least directly, so that it is important to make sure that 
they do not mess up the physics. 

In the modified code we distinguish baryons and CDM 
particles by tagging all of them. We consider the situa- 
tion where 20% of all matter particles are baryonic and 
80% are CDM. At the beginning of each simulation, we 
loop over all particles and for each particle we generate 



a random number from a uniform distribution in [0, 1]. If 
this random number is less than 0.2 then we tag the par- 
ticle as baryon, and otherwise we tag it as CDM. Once 
these tags have been set up they will never been changed 
again, and the code then determines whether the particle 
contributes to the scalar field evolution and feels the fifth 
force or not according to its tag. 

3. Initial Condition and Simulation Parameters 

All the simulations are started at the redshift z = 49. 
In principle, modified initial conditions (initial displace- 
ments and velocities of particles which is obtained given 
a linear matter power spectrum) need to be generated 
for the coupled scalar field model, because the Zel'dovich 
approximation [36, 37] is also affected by the scalar field 
coupling [22]. In practice, however, we have found in our 
linear perturbation calculation [23] that the effect on the 
linear matter power spectrum is negligible (< O(10 -4 )) 
for our choices of parameters 7,/U. Another way to see 
that the scalar field has really negligible effects on the 
matter power spectrum at early times is to look at Fig. 2 
below, which shows that at those times the fifth force is 
just much weaker than gravity and therefore its impact 
ignorable. Considering these, we simply use the ACDM 
initial displacements/ velocities for the particles in these 
simulations, which are generated using GRAFIC2 [38]. 

The physical parameters we use in the simulations are 
as follows: the present-day dark energy fractional energy 
density ^de = 0.76 and Q m = ficDM + = 0.24, Ho = 
73 km/s/Mpc, n s = 0.958, <jg = 0.8. The simulation box 
has a size of 64/i _1 Mpc, where h = i?o/(100 km/s/Mpc). 
We simulate 4 models, with parameters (7, y) equal to 
(0.5, 10" 6 ), (0.5, 10~ 5 ), (1.0, 10" 6 ) and (1.0, 10" 5 ) re- 
spectively (such parameters are chosen so that the de- 
viation from ACDM will be neither too small to be dis- 
tinguishable or too large to be realistic). For each model 
we make 5 runs with exactly the same initial condition, 
but different seeds in generating the random number to 
tag baryons and CDM particles; all the 4 models use the 
same 5 seeds so that results can be directly compared. We 
hope the average of the results from 5 runs could reduce 
the scatter. In all those simulations the mass resolution 
is 1.04 x 10 9 /i -1 Mq, the particle number is 256 3 , the 
domain grid is a 128 x 128 x 128 cubic and the finest 
refined grids have 16384 cells on one side, corresponding 
to a force resolution of ~ 12/i _1 kpc. 

We also make a run for the ACDM model using the 
same parameters (except for (j,, 7, which are not needed 
now) and initial condition. 

B. Preliminary Results 

In Table I we have listed some of the main results for 
the 20 runs we have made, from which we could obtain 
some rough idea how the motions of baryons and CDM 
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FIG. 1: (Colour Online) The snapshots of the distribution of particles. The 4 columns are for the 4 models we consider and the 
3 rows are for three output times, respectively a = 0.2, 0.5 and 1 (corresponding to redshifts 4, 1 and 0) where a is the cosmic 
scale factor. The red crosses denote baryons while blue dots represent CDM particles. All particles are take from a slice of the 
simulation box with 32.0h~ 1 Mpc < z < 32.1/i~ 1 Mpc and projected to an x — y plane. 



particles differ from each other. We see that for the model 
7 = 1.0, \i = 10~ 5 the CDM particles could be up to - 1.6 
times faster than baryons, thanks to the enhancement 
by the fifth force. We will come back to this point later 
when we argue for the necessity of a modified strategy of 
identifying virialized halos. 



TABLE I: The mean velocity of all particles («), baryons («b) 
and CDM particles (wcdm) at z = for the 20 coupled scalar 
field runs. The unit is proper km/s. 





parameters 




run no. 


V 


VB 


VCDM 










1 


367.68 


337.71 


375.18 










2 


367.70 


337.69 


375.19 


7 


= 0.5, /i = 


10 


-6 


3 


367.75 


337.63 


375.28 










4 


367.67 


337.76 


375.15 










5 


367.65 


337.59 


375.16 










1 


382.12 


335.91 


393.68 










2 


382.01 


335.79 


393.55 


7 


= 1.0, /i = 


10 


-6 


3 


381.91 


335.79 


393.55 










4 


382.01 


335.87 


393.54 










5 


382.05 


335.83 


393.61 










1 


412.75 


357.36 


426.60 










2 


412.75 


357.49 


426.55 


7 


= 0.5, n = 


10 


-5 


3 


426.60 


357.44 


438.03 










4 


412.92 


357.62 


426.75 










5 


412.73 


357.44 


426.56 










1 


565.22 


388.77 


609.35 










2 


565.13 


388.68 


609.19 


7 


= 1.0, fj, = 


10 


-5 


3 


565.19 


388.82 


609.30 










4 


565.58 


388.25 


609.66 










5 


564.64 


388.75 


608.63 
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In Fig. 1 we have shown some snapshots of the distri- 
bution of baryonic and CDM particles, to give some idea 
about the hierachical structure formation and for com- 
parisons with other figures below. It shows clearly how 
some clustering objects develop, with filaments connect- 
ing them together. The baryons roughly follow the clus- 
tering of CDM particles, but in some low density regions 
they become slightly separated. 

To understand how the motion of the particles is al- 
tered by the coupling to the scalar field, in Fig. 2 we have 
shown the correlation between the magnitudes of the fifth 
force and gravity on the CDM particles (remember that 
baryons do not feel the fifth force). Ref. [23] has made a 
detailed qualitatively analysis about the general trend of 
this correction, and here we just give a brief description: 

From Eqs. (31, 32) we could see that, when the scalar 
field potential, i.e., the last term of Eqs. (31, 82) could be 
neglected, then the scalar field tp is simply proportional to 
the gravitational potential $ and as a result Eq. (30) tells 
us that the strength of the fifth force is just 2^ 2 times that 
of gravity; in other words, the effect of the scalar field is a 
rcscaling of the gravitational constant by 1 4- 2-f 2 . This is 
because in this situation, the effective mass of the scalar 
field, which is given by rri 2 ^ = d 2 V e f / / dip 2 , where 

V eff (tp) = V{<p) + C%)p DM (38) 

is the effective total potential, is light and the fifth force 
is long-range, like gravity. For comparison, in Fig. 2 we 
also plot this 2j 2 proportion between the two forces, as a 
straight line: lg(/) = lg(0.8g) + lg(27 2 ), where /, g denote 
respectively the magnitudes of the fifth force and gravity, 
and the factor 0.8 in front of g comes from the fact that 
only 80% of the particles are CDM (and thus contribute 
to the fifth force). This scaling relation actually sets an 
upper limit on how strong the fifth force could be relative 
to gravity, should it not be suppressed by other effects. 

In contrast, when the value of ip is small, the last term 
of Eqs. (31, 32) is not negligible and the scalar field ac- 
quires a heavy mass, making it short-ranged. As a result, 
a particle outside a high density region might not feel the 
fifth force exerted by particles in that region, even it is 
quite close to the region. But because it can feel gravity 
from that region, so the total fifth force on the particle 
becomes less than the 2j 2 scaling. 

In general, the value of tp is determined by fj,, 7, pcdm 
as well as its background value (p (which sets the bound- 
ary condition to solve the interior value). At early times 
tp is very close to and pcdm is high everywhere, making 
tp small everywhere too, and suppressing the fifth force so 
that it is significantly below the 27 2 -scaling (first row of 
Fig. 2). At later times, tp increases and pcdm decreases, 
weakening the above effect so that the fifth force becomes 
"saturated" (i.e., approaches the 27 2 prediction) and the 
points in the figure hit the straight lines (last two rows). 
Because decreasing /i and increasing 7 have the same ef- 
fects of making tp small, in the models with /i = 10~ 6 
the fifth force saturates later than in the models with 
/.t = 10~ 5 . In addition, because high pcdm tends to de- 



crease tp and increase the scalar field mass, so in high 
density regions (where gravity is stronger) the fifth force 
also saturates later. 

The agreement between the numerical solution of the 
fifth force and the 27 2 -scaling relation in cases of weak- 
chameleon effect serves as an independent check of our 
numerical code. 

Fig. 3 plots the spatial configuration for the gravita- 
tional potential at the same position and output times 
as in Fig. 1. As expected, the potential is significantly 
deeper where there is significant clustering of matter 
[cf. Fig. 1]. 

We also show in Fig. 4 the spatial configuration for the 
scalar field ip at the same output position and times. At 
early times when ip is small and the scalar field mass is 
heavy, the fifth force is so short-ranged that tp only de- 
pends on the local density. This means that the spatial 
configuration of tp in this situation could well reflect the 
underlying dark matter distribution, a fact which could 
be seen clearly in the first row. As time passes by, the 
mass of the scalar field decreases on average and tp in- 
creases, the value of tp at one point is more and more 
influenced by the matter distribution in neighboring re- 
gions, and such an "averaging" effect weakens the con- 
trast and makes the plots blurring (last two rows). Fur- 
thermore, obviously decreasing [i and increasing 7 could 
increase the scalar field's mass, shorten the range of the 
fifth force, make tp less dependent on its value in neigh- 
boring regions, and thus strengthen the contrast in the 
figures. 



C. Power Spectrum 

To have a more quantitative description about how the 
matter clustering property is modified with respect to 
the ACDM model, we consider the matter power spectra 
P(k) in our simulation boxes. 

The nonlinear matter power spectrum in the present 
work is measured using POWMES [40], which is a public 
available code based on the Taylor expansion of trigono- 
metric functions and yields Fourier modes from a number 
of fast Fourier transforms controlled by the order of the 
expansion. We also average the results from the 5 runs 
for each model and calculate the variance. 

In Fig. 5 shown are the fractional differences between 
the P(k) for our 4 models and for ACDM, at two different 
output times. At early times (a = 0.5, upper solid curves 
in each panel), the difference is generally small, but still 
the 2 models with /1 = 10~ 5 show up to ~ 25% deviation 
from ACDM prediction. This is because for larger fj, the 
scalar field is lighter and the fifth force less suppressed, its 
influence in the structure formation therefore enhanced. 
Notice that on small scales the deviation from ACDM de- 
creases, which is a desirable property of chameleon mod- 
els which are designed to suppress the fifth force on small 
scale high density regions. 

The lower solid curves in each panel of Fig. 5 display 
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FIG. 2: The correlation of the magnitudes of the fifth force and gravity on the CDM particles. The dots denote CDM particles 
taken from a slice of simulation box with 32.0h~ 1 Mpc < z < 32.1/i _1 Mpc, and the lines are the predicted relation between 
fifth force and gravity should the former be not suppressed by the chameleon effect. The 4 columns are for the 4 models we 
consider and the 3 rows are for three output times, respectively a = 0.2, 0.5 and 1 (corresponding to redshifts 4, 1 and 0) where 
a is the cosmic scale factor. 



the same quantities at a = 1.0 (late times). We can see 
the trend of increasing deviation from ACDM for all 4 
models, because fifth force is essentially unsupprcssed at 
the late epoch [cf. Fig. 2]. For example, the deviation of 
the model 7 = 1.0, /i = 10~ 5 is significantly larger than 
that of the model 7 = 0.5, /1 = 10 -5 as naively expected, 
thanks to the lack of suppression of fifth force in both 
models (the 7 = 0.5 model obviously has a smaller satu- 
rated fifth force). 

For comparison wc also plot the AP/P that is pre- 
dicted by the linear perturbation theory for the 4 mod- 
els under consideration (the dashed curves). As can be 
seen there, at large scales (small k) where linear pertur- 
bation is considered as a good approximation, the lin- 
ear and nonlinear results agree pretty well (especially 
for the a = 0.5 case). The largest scale we can probe 
is limited by the size of our simulation boxes (64 Mpc/h) 
and as a result we cannot make plot beyond the point 
k ~ 0.1 h/Mpc, where nonlinearity is expected to first 
become significant. However, in the case of a = 1.0, we 
can see the clear trend of the linear and nonlinear results 
merging towards k < 0.1 h/Mpc at vanishing AP/P. 
Similar results can be found in Fig. 2 of [26] for f(R) 
gravity. 

Because we have two species of matter particles, one 



uncoupled to the scalar field, we are also interested in 
their respect power spectrum and the bias between them. 
These are displayed in Fig. 6. The results could be under- 
stood easily: because a CDM particle always feel stronger 
total force than a baryon at the same position, so the 
clustering of the former is identically stronger than the 
latter as well. This could result in a significant bias be- 
tween these two species at the present time, especially for 
the models with /j, = 10~ 5 , where the fifth force is less 
suppressed. 



D. Mass Functions 

We identify halos in our iV-body simulations using 
MHF (MLAPM Halo Finder) [41], which is the default 
halo finder for MLAPM. MHF optimally utilizes the re- 
finement structure of the simulation grids to pin down 
the regions where potential halos reside and organize 
the refinement hierarchy into a tree structure. Because 
MLAPM refines grids according to the particle density 
on them, so the boundaries of the refinements are sim- 
ply isodensity contours. MHF collect the particles within 
these isodensity contours (as well as some particles out- 
side). It then performs the following operations: (i) as- 
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FIG. 3: (Colour Online) The gravitational potential on the z plane where z = 32h~ 1 Mpc. Dark regions are where the potential 
is more negative (deeper) while light regions are where it is less negative (shallower). The 4 columns are for the 4 models we 
consider and the 3 rows are for three output times, respectively a = 0.2, 0.5 and 1 (corresponding to redshifts 4, 1 and 0) where 
a is the cosmic scale factor. 



suming spherical symmetry of the halo, calculate the es- 
cape velocity v esc at the position of each particle, (ii) if 
the velocity of the particle exceeds v esc then it does not 
belong to the virialized halo and is removed, (i) and (ii) 
are then iterated until all unbound particles are removed 
from the halo or the number of particles in the halo falls 
below a pre-defined threshold, which is 20 in our simula- 
tions. Note that the removal of unbound particles is not 
used in some halo finders using the spherical overdensity 
(SO) algorithm, which includes the particles in the halo 
as long as they are within the radius of a virial density 
contrast. Another advantage of MHF is that it does not 
require a pre-defined linking length in finding halos, such 
as the friend-of-friend procedure. 

When it comes to our coupled scalar field model, the 
MHF algorithm also needs to be modified. The reason is 
that, as we mentioned above, the scalar field ip behaves as 
an extra " potential" (which produces the fifth force) and 
so the CDM particles experience a deeper total "gravita- 
tional" potential than what they do in the ACDM model. 
Consequently, the escape velocity for CDM particles in- 
creases compared with the Newtonian prediction. This is 
indeed important to bear in mind because, as we have 
seen in Table I, the CDM particles arc typically much 
faster than what they are in the ACDM simulation, and 
so if we underestimate v esc then some particles which 



should have remained in the halo would be incorrectly 
removed by MHF. In general, similar things should be 
taken care of in other theories involving modifications to 
gravity in the non-relativistic limit, such as MOND and 
f(R) gravity. 

An exact calculation of the escape velocity in the cou- 
pled scalar field model is obviously difficult due to the 
complicated behaviour of the scalar field, and thus here 
we introduce an approximated algorithm, which is based 
on the MHF default method [42], to estimate it. 

MHF works out v esc using the Newtonian result v\ sc = 
2|$| in which $ is the gravitational potential. Under the 
assumption of spherical symmetry, the Poisson equation 
V 2 <I> = AnGpm could be integrated once to give 

f = (39 ) 
dr r l 

which is just the Newtonian force law. This equation can 
be integrated once again to obtain 

$(r) = G f M( f/V + $ (40) 
Jo r 

where <£>o is an integration constant and can be fixed [42] 
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FIG. 4: (Colour Online) The value of the scalar field (p on the z plane where z = 32/i _1 Mpc. Plotted is \og{yJHip) for better 
contrast. Dark regions are where ip is closer to while light regions are where it is closer to the background value. The 4 columns 
are for the 4 models we consider and the 3 rows are for three output times, respectively a — 0.2, 0.5 and 1 (corresponding to 
redshifts 4, 1 and 0) where a is the cosmic scale factor. 



by requiring that <I>(r = oo) = as 



$n = 



GM V 

Rvii 



G 



R v 



M(< r') 



dr', (41) 



in which R v i r is the virial radius of the halo and Af„j r is 
the mass enclosed in R V i r . 

When the fifth force acts on CDM particles, the force 
law Eq. (39) is modified and these particles feel a larger 
total "gravitational" potential. To take this into account, 
we need to have the knowledge about how the force law is 
modified and a simple rescaling of gravitational constant 
G has been shown to be not physical in certain regimes. 

To solve the problem, we notice that in the MHF code 
Eq. (39) is used in the numerical integrations to obtain 
both $(r) and <J> [cf. Eqs.(40, 41)]. More explicitly, the 
code loops over all particles in the halo in ascending or- 
der of the distance from the halo centre, and when a 
particle is encountered its mass is uniformly distributed 
into the spherical shell between the particle and its pre- 
vious particle (the thickness of the shell is then the dr of 
the integration). Obviously when the fifth force is added 
into Eq. (39) we could use the same method to com- 
pute the total "gravitational" potential which includes 
the contribution from ip (call this contribution be- 
cause dQ^/dr = fifth force, so from Eq. (26) we can eas- 



ily see $ v = -ff-aip). But now there is a subtlety here: 
not all particles arc CDM while only CDM particles con- 
tribute to So in the modified MHF code we calculate 
$,$o and $^,$^0 separately, using all particles for the 
former and only CDM for the latter. Finally 



2|$ + $„-$ -* 



¥>0| 



(42) 



is our estimate of the escape velocity. Because we have 
recorded the components of gravity and the fifth force 
for each particle in the simulation, so the fifth-force-to- 
gravity ratio can be computed at the position of each par- 
ticle, which is approximated to be ($ v — & v o) I — 3>o) 
at that position. In this way we have, at least approxi- 
mately, taken into account the environment dependence 
of the fifth force and thus of <& v . 

To see how our modification of the MHF code affects 
the final result on the mass function, in Fig. 7 we com- 
pare the mass functions for the model 7 = 1.0, /i = 10~ 5 
calculated using three methods: the modified MHF, the 
original MHF (by default for ACDM simulations) and an- 
other modified MHF in which we set v esc to be very large 
so that no particles will ever escape from any halo (this is 
similar to the spherical overdensity algorithm mentioned 
above) . It could be seen that the difference between these 
methods can be up to a few percent for small halos where 
the potential is shallow, particle number is small and the 
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FIG. 5: The fractional difference between the nonlinear matter power spectra for the inhomogeneously coupled scalar field 
model and the ACDM paradigm, where AP(k) is their difference. The range of k is 0.1 to 50 h Mpc -1 . The solid curves are 
mean and the light-blue error bars the standard deviation calculated from the 5 samples of measurements. Upper solid curves 
in each panel: the result when a = 1.0 (redshift 0); Lower solid curves in each panel: the result for a = 0.5 (redshift 1). The 
dashed curves are the corresponding results from the linear perturbation theory for comparison. 



result is sensitive to whether more particles are removed. 
As expected, the second method gives least halos because 
there v esc is smallest and many particles are removed, 
while the third method gives most halos because no par- 
ticles are removed at all. 

Fig. 8 displays the mass functions for the 4 models 
as compared to the ACDM result. As expected, the fifth 
force enhances the structure formation and thus produces 
more halos in the simulation box. 

E. Halo Profiles 

As our simulations have reached a force resolution of 
~ 12/i _1 kpc, while the typical size of the large halos in 
the simulations is of order Mpc, we could ask what the 
internal profiles of the halos look like and how they have 
been modified by the coupling between CDM particles 



and the scalar field. 

We have selected 2 typical halos from each simulation. 
Halo I is centred on (x, y, z) - (27.0, 32.2, lT.S)^ 1 Mpc, 
which is slightly different for different simulations, and 
has a virial mass M v r lr ~ 1.16 x 10 14 /i _1 Mq; halo II is 
centred on (x,y,z) ~ (52.5, 62.9, 38. 3)/i _1 Mpc, which is 
also slightly different for different simulations, and has a 
virial mass M vir ~ 1.70 x 10 13 /i _1 Mq (note here that 
the virial masses are for the A CDM simulations, and for 
scalar field simulations they can be, and generally are, 
slightly different). 

The largest halos, such as halo I, generally reside in the 
higher density regions where the scalar field has a heavier 
mass and shows stronger chameleon effect. Consequently 
the fifth force inside them is severely suppressed so that 
we can expect small deviation from the ACDM halo pro- 
file. On the other hand, the intermediate and small halos 
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FIG. 6: The nonlinear power spectra for CDM and baryonic particles. Shown are the results at two different output times: 
a — 0.5 (dashed curves) and a = 1.0 (solid curves). The upper curves are for CDM particles the lower curves for baryons. Note 
that we only plot the results in the range 0.1 h Mpc -1 < k < 10 h Mpc~ due to the small number of baryonic particles. Note 
also that all the curves are the mean of 5 relevant runs, but we have not included the standard deviations as in other plots, 
because they are small and only make the curves unrecognizable. 



(such as halo II) are mostly in relatively low density re- 
gions in which scalar field has a light mass and the fifth 
force tends to saturate; this means that they should gen- 
erally have a higher internal density than the same halos 
in the ACDM simulation due to the enhanced growth by 
the fifth force. 

This above analysis is qualitatively confirmed by Fig. 9, 
where we can see that for the cases of fj, = 10~ 6 (stronger 
chameleon) the difference between the predictions of the 
coupled scalar field models and the ACDM paradigm is 
really small for halo I. Furthermore, this figure also shows 
some new and more interesting features for the cases of 
fj, = 10 -5 . Considering halo I in the models 7 = 0.5, /1 = 
10~ 5 and 7 = 1.0, /i = 10~ 5 for example: in the former 
case, the coupled scalar field model produces an obviously 
consistent higher internal density profile than ACDM, 
from the inner to the outer regions of the halo; while for 
the latter case, the density profile of the coupled scalar 



field model is lower in the inner region but higher in the 
outer region! 

We plan to make a detailed analysis of the complexities 
arising here regarding the effects of a coupled scalar field 
on the internal density profiles for halos in a future work, 
and in this work we will only give a brief explanation for 
the new feature observed above: Fig. 10 shows the dis- 
tributions of particle velocities and fifth-force-to-gravity 
ratio at the positions of the particles in halo I for the two 
models 7 = 0.5, /1 = 10~ 5 and 7 = 1.0,^ = 10 -5 . As we 
can see there, for the model 7 = 1.0, fx = 10~ 5 the large 
value of 7 makes the chameleon effect strong so that the 
fifth force is generally much smaller than gravity in mag- 
nitude, while at the same time the velocities of particles 
are more concentrated towards the high end, implying a 
significantly higher mean speed than ACDM (as we have 
checked numerically) ; as a result in the central region of 
the halo the particles have higher kinetic energy than in 
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FIG. 7: The comparison of three different methods to identify halos from the simulations. The first method (solid curve) uses 
the approximation Eq. (42) to estimate the escape velocity; the second method (dotted curve) uses the ACDM prediction, 
which underestimates v eac and causes more particles to be removed; the third method (dashed curve) uses v esc — oo so that 
no particles are removed. All curves are for the model 7 = 1.0, [i — 10~ J , and the curves and error bars denote respectively the 
mean and standard deviation for 5 sample results. 



ACDM but the potential is not significantly deeper, so 
that particles tend to get far away from the centre of the 
halo, producing a lower inner density profile. As for the 
model with 7 = 0.5, /1 = 1CP 5 , the situation is just the 
opposite: the fifth force is saturated and of the same order 
as gravity due to the weak chameleon effect, so that the 
total potential is greatly deeper than its ACDM counter- 
part, while at the same time the particles are not as fast 
as those in 7 = 1.0, fi = 10 -5 : the consequence is that 
the halo accretes more particles towards its centre. 

One might also have interests in how the CDM part 
and the baryonic part of the halo profile differ from each 
other, and this is given in Fig. 11, in which we compare 
the CDM and baryon density profiles of halos I and II. 
Obviously the CDM density is higher than baryons ev- 
erywhere, again thanks to the boost from the fifth force. 
For smaller fi and larger 7 (which help strengthen the 
chameleon effect and suppress the fifth force) , the differ- 
ence between the two is smaller (compare the halo I in 
the models 7 = 0.5, fi = 10~ 5 and 7 = 1.0, /1 = 10~ 5 ). 
However, in the situations where the fifth force is already 
saturated or unsuppressed (such as in halo II) , increasing 
7 increases the saturated value of the fifth force, and thus 
can instead magnify the difference (compare the halo II 
in the models 7 = 0.5, fi = 10~ 5 and 7 = 1.0, [i = 10 -5 ). 



The above results again show the complexity of the 
chameleon scalar field model as compared to other cou- 
pled scalar field models. We will study the environment 
and epoch dependence of the halo density profiles in a 
upcoming work. 



IV. DISCUSSION AND CONCLUSION 

To summarize, in this paper we have investigated into 
a model where CDM and baryons couple differently to 
a chameleon-like scalar field, and performed full A^-body 
simulations, by directly solving the spatial distribution 
of the scalar field, to study the nonlinear structure for- 
mation under this setup. The new complexity introduced 
here compared to [23] is that we must distinguish between 
baryons and CDM so that we know how to calculate the 
force upon each particle. We do this by tagging the initial 
almost homogeneously distributed particle randomly so 
that 80% of all particles are tagged as CDM. We then 
only use CDM particles to calculate the scalar field and 
only applies the fifth scalar force on CDM particles. 

The coupling function (characterized by the coupling 
strength 7) and bare potential (characterized by the pa- 
rameter fi which controls its steepness) of the scalar field 
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FIG. 8: The mass functions for the 4 models under consideration at a = 1.0 (solid curves). The dashed curve is the mass 
function for the ACDM simulation. The (solid) curves and error bars denote respectively the mean and standard deviation for 
5 sample results. 



are chosen to be the same as in [23]. As discussed there, 
the coupling of the scalar field to CDM particles acts on 
the latter a fifth force. When fj, is small and 7 is large, 
the chameleon effect becomes stronger, which gives the 
scalar field a heavy mass, making the fifth force short- 
ranged and the scalar field dependent only on the local 
matter density. Other ways to strengthen the chameleon 
effect includes increasing the density and decreasing the 
background value of the scalar field, which itself is equiv- 
alent to increasing the background CDM energy density. 
We have displayed in Figs. 2, 4 how changing the deter- 
mining factors change the scalar field configuration and 
the strength of the fifth force. 

We have also measured the nonlinear matter power 
spectrum from the simulation results and compared them 
with the ACDM prediction. Depending on the values of 
7, /i as well as the background CDM density, the former 
can be up to ~ 80% larger than the latter. Nonetheless, 
when the chameleon effect is set to be strong, the devia- 



tion gets suppressed and in particular decreases towards 
small scales, showing the desirable property of chameleon 
models that they evade constraints on small scales. The 
bias between CDM and baryons power spectra follows 
the same trend. 

To identify virialized halos from the simulations, we 
have modified MHF, MLAPM's default halo finder, so 
that the calculation of the escape velocity includes the 
effect of the scalar field. We find that such a modifica- 
tion leads to up to a-few-percent enhancement on the 
mass function compared with what is obtained using the 
default MHF code, because in the latter case the escape 
velocity is underestimated and some particles are incor- 
rectly removed from the virialized halos. We find that the 
mass function in the coupled scalar field models is sig- 
nificantly larger than the ACDM result, because of the 
enhanced structure growth induced by the fifth force. 

Finally, we have analyzed the internal profiles of (the 
same) two halos selected from each simulation. We find 
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FIG. 9: The overdensity inside two typical halos selected from the simulations for each model (see text for a more detailed 
description), as a function of the halo radius R. The solid curves and error bars denote respectively the mean and standard 
deviation for the results from our 5 sample runs, while the dashed curve is the ACDM prediction. In each panel the upper two 
curves (black) at large radii are for halo I and the lower two curves (green) at large radii are for halo II. pb is the background 
matter (baryons and CDM) density. We skipped the result for R < 10/i -1 kpc, which is not reliable due to the limit from 
resolution. 



that when chameleon effect is strong and fifth force is 
suppressed, our result is very close to the ACDM predic- 
tion; this is the case for very large halos (which generally 
reside in higher density regions) and n = 10~ 6 . For large 
halos and fi = 10~ 5 , the situation is more complicated 
because the competition between two effects of the cou- 
pled scalar field, namely the speedup of the particles and 
the deepening of the total attractive potential, has ar- 
rived a critical point: if the former wins, such as in the 
model 7 = 1.0, fi = 10~ 5 , then the inner density profile 
can be lower than that in ACDM; if the latter wins, such 
as in the model 7 = 1.0, fx = 10~ 5 , then we expect the 
opposite. For smaller halos which locate in lower den- 
sity regions, the fifth force is not suppressed that much 
and causes faster growth of the structure, so the halo 
is more concentrated and has a higher internal density. 
Meanwhile, the bias between baryons and CDM density 



profiles also increase as the fifth force becomes less sup- 
pressed, which is as expected. 



Our results already show that new features can be 
quantitatively studied with the iV-body method and 
the improving supcrcomputing techniques, and that the 
chameleon model has rather different consequences from 
other coupled scalar field models. The enhancement in 
the structure formation due to the fifth force is signifi- 
cant for some of our parameter space, which means other 
observables, such as weak lensing, could place new con- 
straints on the model. Also one might be interested in 
how the halo profiles would be like at different epoch of 
the cosmological evolution and in different environments. 
These will be left as future work. 
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particles, for the two models 7 = 0.5, fi = 10~ 5 and 7 = 1.0, fi = 10 -5 . Each dots represent a particle. 
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Appendix A: Discretized Equations for Our A-body 
Simulations 



In the MLAPM code the partial differential equation 
Eq. (31) is (and in our modified code Eq. (32) will also 
be) solved on discretized grid points, and as such we must 
develop the discretized versions of Eqs. (29 - 32) to be 
implemented into the code. 

But before going on to the discretization, we need to 
address a technical issue. As the potential is highly non- 
linear, in the high density regime the value of the scalar 
field y/Rtp will be very close to 0, and this is potentially 
a disaster as during the numerical solution process the 
value of y/Hip might easily go into the forbidden region 
ip < [25]. One way of solving this problem is to define 
X = X e " m which x is the background value of X; as hi 
[25]. Then the new variable u takes value in (—00,00) so 
that e u is positive definite which ensures that \ > 0. 
However, since there are already exponentials of x m 
the potential, this substitution will result terms involving 
cxp [cxp(u)], which could potentially magnify any numer- 
ical error in u. 

Instead, we can define a new variable u according to 



1 



(Al) 



By this, u still takes value in (—00, 00), e u <E (0, 00) and 
thus e x € (1, 00) which ensures that x is positive definite 
in numerical solutions. Besides, e^ x = [1 + e v ] so that 
there will be no exponential-of-exponential terms, and 
the only exponential is what we have for the potential 
itself. j3 = — 1 above. 

Then the Poisson equation becomes 
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where we have defined ily = k V{<p)/3Hq which is deter- 
mined by background cosmology, the quantity e 



also determined solely by background cosmology. These 
background quantities should not bother us here. 
The scalar field EOM becomes 
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in which we have used the fact that x — log(l + e u ) =>■ 
= !^ e u an d moved all terms depending only on 
background cosmology (the source terms) to the right 
hand side. 

So, in terms of the new variable u. the set of equations 
used in the TV-body code should be 



dx c 
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(A5) 



plus Eqs. (A2, A3). These equations will ultimately be 
used in the code. Among them, Eqs. (A2, A5) will use the 
value of u while Eq. (A3) solves for u. In order that these 
equations can be integrated into MLAPM, we need to 
discretize Eq. (A3) for the application of Newton-Gauss- 
Seidel iterations. 

To discretize Eq. (A3), let us define b = jj^i- The 
discretization involves writing down a discretion version 
of this equation on a uniform grid with grid spacing h. 
Suppose we require second order precision as is in the 
standard Poisson solver of MLAPM, then Vu in one di- 
mension can be written as 



Vit V n Uj = 



2h 
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where a subscript j means that the quantity is evaluated 
on the j-th point. Of course the generalization to three 
dimensions is straightforward. 

The factor b in V • (bVu) makes this a standard variable 
coefficient problem. We need also discretize b, and do it 
in this way (again for one dimension): 
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where we have defined bj + i = (bj + &j+i) /2 and bj_i = 
(bj-i + bj) /2. This can be easily generalize to three di- 



mensions as 
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Then the discrete version of Eq. (A3) is 
L h (u iijlk ) = 0, 
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Then the Newton-Gauss-Seidel iteration says that we can 
obtain a new (and often more accurate) solution of u, 
u^Yki using our knowledge about the old (and less accu- 
rate) solution uff k as 



The old solution will be replaced by the new solution to 
Uij jk once the new solution is ready, using the red-black 
Gauss-Scidel sweeping scheme. Note that 
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In principle, if we start from a high redshift, then the 
initial guess of Uij & could be such that the initial value 
of x m an the space is equal to the background value x, 
because anyway at this time we expect this to be approx- 
imately true. For subsequent time steps we could use the 
solution for Ui,j t k at the previous time step as our initial 
guess; if the time step is small enough then we do not ex- 



pect u to change significantly between consecutive times 
so that such a guess will be good enough for the iteration 
to converge fast. 

In practice, however, due to specific features and al- 
gorithm of the MLAPM code [32], the above procedure 
may be slightly different in details. 



